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ABSTRACT 

Context. Ionizing radiation plays a crucial role in star formation at all epochs. In contemporary star formation, ionization abruptly 
raises the pressure by more than three orders of magnitude; the temperature increases from ~ 10 K to ~ lO'' K, and the mean molecular 
weight decreases by a factor of more than 3. This may result in positive feedback (either by compressing pre-existing clouds and 
rendering them unstable or by sweeping up gravitationally unstable shells). It may also result in negative feedback (by dispersing 
residual dense gas). Ionizing radiation from OB stars is also routinely invoked as a means of injecting kinetic energy into the interstellar 
medium and as a driver of sequential self-propagating star formation in galaxies. 

Aims. We describe a new algorithm for including the dynamical effects of ionizing radiation in SPH simulations, and we present 
several examples of how the algorithm can be applied to problems in star formation. 

Methods. We use the HEALPix software to tessellate the sky and to solve the equation of ionization equilibrium along a ray towards 
each of the resulting tesserae. We exploit the hierarchical nature of HEALPix to make the algorithm adaptive, so that fine angular 
resolution is invoked only where it is needed, and the computational cost is kept low. 

Results. We present simulations of (i) the spherically symmetric expansion of an Hii region inside a uniform-density, non-self- 
gravitating cloud; (ii) the spherically S5mimetric expansion of an Hii region inside a uniform-density, self-gravitating cloud; (iii) 
the expansion of an off-centre Hn region inside a uniform-density, non-self-gravitating cloud, resulting in rocket acceleration and 
dispersal of the cloud; and (iv) radiatively driven compression and ablation of a core overrun by an Hn region. 
Conclusions. The new algorithm provides the means to explore and evaluate the role of ionizing radiation in regulating the efficiency 
and statistics of star formation. 
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1. Introduction 

The influence of the UV radiation emitted by a massive star (or 
small group of massive stars) on the ambient interstellar medium 
is of great interest. As this radiation ionizes the surrounding neu- 
tral material, it creates a region of hot ionized gas (Str0mgren 
1939), known as an Hn region. The huge contrast in pressure 
between the ionized gas in the Hn region and the surrounding 
neutral gas causes the Hii region to expand and sweep up the 
neutral gas into a dense shell (Kahn 1954; Spitzer 1978). This 
dense shell is likely to be prone to many diflFerent instabilities 
(e.g. Elmegreen 1994; Vishniac 1983), and some of these insta- 
bilities may trigger fragmentation and star formation. This mode 
of star formation is usually referred to as the "collect and col- 
lapse" mechanism (Elmegreen & Lada 1977; Whitworth et al. 
1994a,b). Deharveng et al. (2003, 2005, 2008) and Zavagno et 
al. (2007) have observed several instances where this mechanism 
appears to be operating. 

The structure of the intersteUar medium is observed to be 
extremely irregular and to contain many clouds. As an Hn re- 
gion expands, it may overrun and compress nearby clouds, caus- 
ing them to collapse. This mechanism is known as "radiation 
driven compression", and has been simulated by several work- 



ers. Some of these simulations (Sandford et al. 1982; Bertoldi 

1989; Lefloch & Lazarefi^ 1994; Miao et al. 2006; Henney et 
al. 2008) concentrate on the morphology of the resulting bright- 
rimmed clouds. In contrast, Peters et al. (2008) address the possi- 
bihty that the flow of ionized gas off an irradiated cloud might be 
an important source of turbulence; and Kessel-Deynet & Burkert 
(2003), Gritschneder et al. (2008), Dale et al. (2005) and Miao 
et al. (2008) explore the possibility that the collapse of a bright- 
rimmed cloud might sometimes lead to triggered star forma- 
tion. These simulations show that it is possible to reproduce 
the observed morphologies of bright-rimmed clouds (Lefloch 
& Lazarefi^ 1995; Lefloch et al. 1997; Morgan et al. 2008). 
However, the evolution of bright-rimmed clouds and the role of 
expanding Hn regions in triggering star formation are still poorly 
understood. For example, Dale et al. (2007) argue that the main 
effect of an expanding Hn region may simply be to expose stars 
that would have formed anyway. 

A variety of techniques have been used to follow the 
propagation of ionizing radiation in hydrodynamic simulations. 
Kessel-Deynet & Burkert (2000) consider each SPH particle to 
be at the end of an individual ray, and compute the attenua- 
tion of ionizing radiation along this ray in order to determine 
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whether the particle at the end is ionized. Dale et al. (2005, 
2007) have refined this method so as to make it less compute- 
intensive. Alvarez et al. (2006) have invoked the HEALPi)|J al- 
gorithm in order to distribute rays from the ionizing source in a 
regular manner HEALPix (Gorski et al. 2005) is a versatile tree 
structure which distributes pixels over the celestial sphere, to a 
user prescribed level of refinement. At each successive level of 
refinement a pixel from the level above is split into four smaller 
pixels. Each pixel corresponds to the end of a ray and each ray 
at a given level represents an equal solid angle on the celestial 
sphere. Abel & Wandelt (2002) have shown how the HEALPix 
rays can be split adaptively, to increase the level of refinement 
where this is needed. This adaptive ray splitting has been used 
by Abel et al. (2007) and Krumholz et al. (2007). 

In this paper we simulate the expansion of an Hii region us- 
ing the Smoothed Particle Hydrodynamics (SPH) code SEREN 
(Hubber et al., in preparation). SPH (Lucy 1977; Gingold & 
Monaghan 1977) is a Lagrangian numerical method that de- 
scribes a fluid using particles. The SEREN code is described in 
Section 3. The transport of ionizing radiation uses the hierar- 
chical ray structure of HEALPix, with rays emanating from the 
ionizing source into the surrounding medium. We walk along 
each ray with an adaptive step proportional to the local smooth- 
ing length. At each step we calculate the attenuation of ionizing 
radiation by computing the number of recombinations into ex- 
cited states of H°, i.e. we adopt the on-the-spot approximation 
(Osterbrock 1974). Rays are split wherever the separation be- 
tween particles becomes less than the separation between rays. 
This is a robust criterion, which ensures that a ray is split when 
it encounters a dense region, whilst a neighbouring ray passing 
through a more rarefied region is not split. Once the ionization 
front is located along a ray, the propagation of that ray is ter- 
minated, in order to reduce the computational overhead; this is 
particularly important during the early stages of evolution, when 
the Str0mgren sphere only involves a small fraction of the total 
number of SPH particles. 

Although heavy elements play an important role in determin- 
ing the temperature of interstellar gas, we do not consider their 
chemistry here. Instead, we assume that the composition of the 
gas is Z = 0.7 hydrogen and Y - 0.3 helium, by mass; and 
that the helium is everywhere neutral. In the neutral gas we as- 
sume that the temperature is Tn = lOK and the hydrogen is all 
molecular; hence the mean molecular weight is /i„ = 2.35 and 
the isothermal sound speed is c,, = 0.2 km s"' . In the ionized gas 
we assume that the temperature is = 10"* K and the hydro- 
gen is all ionized; hence jji = 0.678 and q = 11 kms '. In the 
transition region between the molecular and ionized regions, we 
impose a linear temperature gradient between these two limiting 
values. 

We explore three cases, (i) In the first case, the star lies at 
the centre of a spherical cloud, and the Hii region is spherically 
symmetric throughout its expansion, (ii) In the second case, the 
star is placed off-centre inside a spherical cloud; here the Hii 
region breaks out of the cloud on one side, and the remainder 
of the cloud is accelerated by the rocket effect, (iii) In the third 
case, the star is located outside the cloud from the outset, and a 
shock is driven into the cloud ahead of the ionization front (i.e. 
radiation driven compression). These cases are presented only as 
illustrative examples of what the code can simulate. Detailed in- 
vestigations of these phenomena will be presented in subsequent 
papers. 



' http://healpix.jpl.nasa.gov 



The paper is organized as follows. In Section |2] we discuss 
briefly the expansion of an Hii region, once the initial Str0mgren 
sphere has formed. In Section |3] we describe in detail how we 
treat the propagation of ionizing radiation. In Section |4] we test 
the algorithm on the three cases described above, and in Section 
|5]we discuss the results and conclude. 

2. The D-type expansion of an Hii region 

2.1. Spitzer solution 

Consider an arbitrary density field, p(r), and suppose that there 
is an ionizing star at the centre of co-ordinates. Assuming ion- 
ization equilibrium, and neglecting the diffuse radiation field, the 
position of the ionization front (IF), in the direction of the unit 
vector e , is given by R,p = R^^e , where R^^ is obtained from 

p^ire)r'dr^—^^I^,,. (1) 
Jr=o 4;rQ'3 

Here, m - m„ IX = 2.4 X 10^24 g jjjg 

mass associated with 
one hydrogen nucleus when account is taken of the contribution 
from helium, Wp is the proton mass, Ni^^., is the rate at which 
the exciting star emits Lyman continuum photons, and is the 
recombination coefficient into excited states only. Eqn. ([T]l as- 
sumes that the material inside the Hii region is fully ionized. 

Thus Eqn. ([TJ determines the radius R^^^ at which all the ion- 
izing photons emitted in the direction e have been used up bal- 
ancing recombinations into excited states. We ignore recombina- 
tions straight into the ground state by invoking the on-the-spot 
approximation (Osterbrock 1974). 

Str0mgren (1939) was the first to show that the transition 
from a state of almost completely ionized material to a state of 
almost completely neutral material occurs in a very short dis- 
tance compared with the dimensions of the Hii region. For ex- 
ample, for a spherically-symmetric Hii region expanding into a 
cloud of uniform density p^, the radius of the ionization front is 
given by the Str0mgren radius 



LyC 


1/3 

^ 0.3 pc 




Pn 




1049 s-1 J 


\ 10-20 gcm-3j 



and the distance over which the degree of ionization changes 
from 90 % to 10 % is given by 

At;^, . . 0.0002 pc( /" \' . (3) 

PnO" \10 ^Ogcm 

Here cr - 7 x lO"'** cm2 is the mean photoionization cross sec- 
tion presented by a hydrogen atom to Lyman continuum photons 
from an OB star. 

Because the squared sound speed in the ionized gas in- 
side the Hii region is ^ 122km2s-2, whereas the squared 
sound speed in the neutral material outside the Hii region is 
cl ^ 0.0352 km2s-2, there is a large pressure difference be- 
tween the two regimes (more than three orders of magnitude), 
and this results in rapid expansion of the Hii region. The out- 
ward propagation of the ionization front is subsonic relative to 
the ionized gas (where the sound speed is c; ~ 11 kms"'), but 
supersonic relative to the neutral gas (where the sound speed 
is Cn ~ 0.2 km s"'). Consequently a strong shock front appears 
ahead the ionization front (for a detailed discussion see Kahn 
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1954). Spitzer (1978) has obtained an approximate analytic so- 
lution for this phase of the evolution. In this solution, the radius 
of the ionization front is given by 



the density of the ionized gas by 



pi = Pn 1 + 



7 Cit 
4^ 



-6/7 



the mass of ioruzed gas by 



Mi(f) = 



1 + 



Q'b Pn \ 4 /Jj, 

and the speed at which the ionization front propagates by 



= c, 1 + 



7 Cit 
4^ 



-3/7 



(4) 



(5) 



(6) 



(7) 



2.2. Semi-analytic approacln 

A similar result can be obtained by solving the equation of mo- 
tion of the accreting thin shell accelerated by the thermal pres- 
sure of the Hn region (see e.q. Hosokawa & Inutsuka 2006). We 
use this approach to obtain a more general solution which in- 
cludes effects of the gravity and use this solution to calibrate our 
numerical code with the self-gravity switched on (see ^4.2| i. 

The equation of motion of the shell of mass Msh and radius 
^sh is 



dt 



{MM = AnR\P; 



GM,y,Mi GMl 



2R\ 

sh 



(8) 



vl/2 



Where P, = p^cf, M, = (4;r/?>i)/3 and = (S^) 

the pressure, the mass and the density of the ionized gas, respec- 
tively. The second term on the rhs represents the gravitational 
force due to the Hii region mass, and the third term is the aver- 
age gravitational force acting on the shell due to its own mass 
(it is zero at its inner edge and GM^JR^^^ at the outer one; see 
Whitworth & Francis 2002). 

The mass of the shell is a difference between the mass of the 
original neutral gas which occupied the sphere of the radius ^sh 
and the mass of the Hii region 



471 

Msh = y^sh(Pn-Pi) • 



(9) 



Inserting (j9]l into (j8]l and using Eqn. (j2]i we get 



In 



-3ci«'-yGp„< 



St '"^sh 

1-1^ 

^sh 



(10) 



We solve this equation numerically using the subroutine 
integrate . odeint from the Scientific python library (Jones et 
al. 2001). Since the model assumes an infinitesimally thin shell, 
the shell radius /?sh can be compared to the average radius of the 
ionization front and the average radius of the shock front /?sf in 
the numerical simulations (see Fig. 1 1 and related discussion in 



3. Numerical treatment 

We use the SPH code SEREN which has been extensively 
tested using a wide range of standard and non-standard prob- 
lems. SEREN is an astrophysical SPH code parallelized us- 
ing OpenMP and designed primarily for star and planet for- 
mation problems. It contains both the traditional SPH method 
(Monaghan 1992) and the more recent 'grad-h' SPH formula- 
tion (Price & Monaghan 2004). There exists the option to in- 
clude additional features within the basic SPH algorithm, such 
as time-dependent and/or Balsara-switched viscosity (Morris & 
Monaghan 1997; Balsara 1995), artificial conductivity (Price 
2008) and a variety of equations of state. The SPH equations 
can be solved with a choice of integrators, such as the symplec- 
tic 2nd order Leapfrog integrator, in conjunction with a block 
time-stepping scheme. 

SEREN uses a Barnes-Hut octal-spatial tree (Barnes & 
Hut 1986), both to speed up the calculation of gravitational 
forces, and to facilitate the collation of particle neighbour lists. 
The Barnes-Hut tree calculates the gravitational force up to 
quadrupole or octupole order, and has the option of using ei- 
ther the traditional geometrical opening angle criterion or the 
GADGET-style multipole acceptance criterion (Springel et al. 
2001). For problems where multiple gravitationally bound ob- 
jects form, SEREN uses sink particles (Bate et al. 1995) to fol- 
low the simulation beyond the formation of the first dense bound 
objects. A 4th-order Hermite integrator (Makino & Aarseth 
1992) is included, so that the code can follow the later N-body 
evolution of the sinks, after the accretion phase has finished. 
A paper describing SEREN in detail, and the tests to which it 
has been subjected, will be submitted shortly (Hubber et al., in 
preparation). 

3.7. Ray casting 

In order to determine the overall shape of the ionization front 
we consider many different directions e . These directions are 
chosen using the HEALPix algorithm (Gorski et al. 2005). 
HEALPix generates a set of A^f - 12 x 4^ directions (hereafter 
'rays'), where { is the level of refinement and ^ = is the lowest 
level. The rays of level f are distributed uniformly over the ce- 
lestial sphere, as seen from the ionizing star Each ray on level i 
intersects the celestial sphere at the centre of an approximately 
square element of solid angle 

An TT 
AQ/ = — steradians = r steradians. 

Nc 3 X 4^ 

Therefore, the angle between neighbouring rays is of order 
A0{ (AQty^ ^ 2"^ radians. 

Wlien HEALPix goes to a higher level, the angular res- 
olution is increased by splitting each of the rays at the previous 
level, {- 1, into four rays (see ^3.2\ . In the simulations presented 
here, we limit the ray splitting to ( - 1, so there is a maximum 
of 196,608 rays at the ionization front. 

The integral on the left-hand side of Eqn. ([T]) must be evalu- 
ated numerically along each ray until it equals 7^,^^ . To perform 
this integration we define a set of discrete evaluation points along 
each ray. At each evaluation point, j, we determine the density pj 
using an SPH summation over the A/j^Em nearest SPH particles, 
k. 



Pj 




(11) 
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The integral in Eqn. ([T]l is then evaluated using the trapezium 
method 



I{r 



Jr=0 

j'=l 



fihf-i 



(12) 
(13) 



where fi is a dimensionless tolerance parameter of order unity 
controlling the integration step, and hf is the smoothing length 
of evaluation point /. The identifier / = denotes the eval- 
uation point placed at the position of the star, and consecutive 
evaluation points are then located according to 



+ flhyt. 



(14) 



Thus fihji represents the adaptive step along the line of sight 
from the ionizing star towards the ionization front. Fig. [T] shows 
a representation of the method. Acceptable accuracy is obtained 
with fx - 0.25 and this is the value we use in the simulations 
presented here. 

3.2. Ray splitting 

The use of ray splitting allows us to maintain good resolution 
while achieving significant speed-up in comparison with uni- 
form ray tracing (cf. Abel & Wandelt 2002). In our scheme a ray 
is split into four child-rays as soon as its linear separation from 
neighbouring rays, rjlS.Q(, exceeds fihj. Here Vj is the distance 
from the ionizing star to evaluation point j, fi is a dimensionless 
parameter controlling the angular resolution of the ensemble of 
rays representing the ionizing radiation, and hj is the smoothing 
length of the local evaluation point /'. Hence we require 



t > logj 



fihj 



(15) 



Acceptable results are obtained with fi in the range (1.0, 1.3). 
A smaller fi gives greater accuracy, but at the expense of more 
HEALPix rays. As one moves outwards away from the central 
star, rays are only ever split, they are never merged. 

The NESTED version of HEALPix allocates identifiers to rays 
according to the scheme illustrated in Fig. |2] which shows a 
small patch on the celestial sphere. Here the mother-ray of ray m 
has ID 



- -int(-) 

and the four child-rays of ray m have IDs 
m" - Am + n, with n = 0, 1 , 2, 3 



(16) 



(17) 



Here n - Q corresponds to the child-ray to the South, n - I 
corresponds to the child-ray to the West, n - 2 corresponds to 
the child-ray to the East and n = 3 corresponds to the child-ray 
to the North. This numbering scheme makes it very simple to 
trace a ray back to the star at the origin, provided you know the 
ID of the ray and its level. 

3.3. Ray rotation 

To avoid numerical artifacts due discretization in the directions 
of rays, we rotate the ensemble of HEALPix rays through three 




Fig. 2. Illustration of the scheme used for splitting a ray using 
HEALPix's NESTED scheme. This map shows the tessellation 
of the celestial sphere corresponding to three successive levels 
of HEALPix. Each tessera has a ray at its centre. The big bold 
tessera represents the solid angle AOf of a single ray at level I. 
The four intermediate tesserae represent the solid angles of its 
child-rays at level £ +\. And the sixteen smallest tesserae repre- 
sent the solid angles of its grandchild-rays at level £ + 2. 



random angles (about the z-, x'-, and z"-axes), each time we 
build the ensemble. This process is necessary to prevent the for- 
mation of artificial corrugations in the ionization front, as ex- 
plained in detail in Krumholz et al. (2007). 



3.4. Setting temperatures 

At each evaluation point we estimate the value of I{rj) using 
Eqn. 



13 I. If the calculation returns I(rj) < Z^^^, then the evalu- 



ation pomt lies in the interior of the Hn region . Therefore, some 
of the particles that belong to its neighbour list should be ion- 
ized. We flag a particle as ionized, if it passes the following two 
tests: (a) its coordinates are inside the solid angle of the ray; and 
(b) its distance from the star is smaller than the distance of the 
evaluation point (rj). If a particle is ionized, it is given a tem- 
perature T - Ti. Otherwise it remains neutral with temperature 
T - Tn. Fig|3]illustrates the location of the neighbours of evalu- 
ation point j, and the ones which pass these tests. 




Fig. 3. This figure shows a single ray (soUd line) along with the 
boundaries of its solid angle (dash lines); the star lies on the left 
hand side and the ionization front on the right. The big black dot 
is the evaluation point j, and the circle is its smoothing region. 
The dot dashed line shows the position where the distance from 
the star is the same as that of the evaluation point. Small grey 
dots are particles which are flagged as being ionized, and the 
crosses are particles which are flagged as being neutral. 
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Fig. 1. This figure shows the evaluation points (big black dots) along a family of rays (solid lines), starting from the ionizing star. 
To avoid confusion, the ray casting is plotted in 2-D. A ray is split as soon as its linear separation from neighbouring rays, rjAOc, 
exceeds fJ^j, where hj is the smoothing length of the local evaluation point. Evaluation points where rays are split have been marked 
with open circles to indicate that they make no contribution to the summation in Eqn. ( 13 i. A binary chop subroutine is used to 
locate the position of the ionization front (see 5 3.5 i. 



3.5. Location of the ionization front 



13 I and 



We store in memory the value of the integral in Eqn. 
the coordinates of the evaluation point at the end of a given ray 
segment. The distance R^^ of the ionization front from the star 
along a given ray is given by the condition /(/?,p) = /^ax- 
is determined using a binary chop algorithm between the last 
two evaluation points ry_i and rj for which /(ry_i) < < 
/(rj).The binary chop algorithm is stopped once the adjustment 
to rj becomes smaller than lO^^/iy-i. 

If the material of the cloud along the line of sight of a ray is 
sufficiently rarefied, / (r) may never reach /^^^ within the com- 
putational domain. We therefore estimate the maximum distance 
Oj.max of all gas particles and if rj > r^ ^ax we stop calculating 
the sum. Any ray that satisfies this condition is characterized as 
open and the position of the ionization front is not determined 
along such rays (see Fig.Hl). 



3.6. Temperature smootiiing 

It is not feasible to resolve the ionization front in a standard 
SPH simulation of an evolving Hn region (see Appendix |A|i. 
Consequently we have to smooth the temperature gradient ar- 
tificially across the ionization front. To demonstrate why this is 
necessary, we model a uniform density spherical cloud consist- 
ing of 10^ pre-settled SPH particles. The total mass of the cloud 
is M = IOOOMq, and the radius is R - Ipc. At the centre 
of the cloud there is a single star emitting A/'lc = 10'*'' ioniz- 




Fig. 4. This figure shows how a ray penetrates the borders of the 
cloud. Once it finds no more neutral material to ionize along its 
line of sight, it is characterized as open. 



ing photons per second. There are no gravitational forces. We 
simulate this system in two ways: in Case A, the temperature 
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Fig. 5. The smoothing length of all (10^) particles are plotted as 
a function of their distance from the star at f = 0.14Myr. Each 
point represents an individual SPH particle. The ionization front 
is located at r,p ~ 0.63 pc away from the star. A gap between the 
ionized and neutral regimes is formed due to the discontinuous 
transition in temperature (Case A). 



Fig. 6. The smoothing lengths of the particles are plotted as a 
function of their distance from the star, for the Case B. Despite 
that the snapshot is taken at the same time as in Fig. [5] (f = 
0. 14 Myr) the position of the ionization front is placed further, at 
~ 0.87 pc. The transition between the two regimes is contin- 
uous now. 



gradient is unsmoothed; and in Case B, the temperature gradient 
is smoothed. 

In Case A, the temperatures of the SPH particles are set as 
described in 5 3.4 and therefore there is an abrupt decrease in 
temperature across the ionization front. The simulation is termi- 
nated once the ionization front reaches the edge of the cloud. 
In Fig. [5] the smoothing lengths of all SPH particles are plot- 
ted against their distance from the ionizing source at time t - 
0.14 Myr. Due to the abrupt temperature and pressure change 
across the ionization front, a gap is formed between the ionized 
and neutral regions. This gap prevents the ionization of new ma- 
terial from the neutral region during the expansion, and therefore 
the mass of the Hii region does not increase (see Fig.]?}. This in 
turn alters the evolution of the radius of the Hii region, as shown 
in Fig. [8] 

In Case B, the temperature gradient across the ionization 
front is smoothed. Well inside the Hn region, and well outside 
the ionization front, we set the temperatures of the SPH parti- 
cles as described in §3.4 However, in a layer of thickness 2h^^ 
across the ionization front we linearly interpolate between the 
temperatures in the ionized and neutral regimes, i.e. 



T(r) = ^ (r„ + Ti) + '—^ (r„ - Ti) , \r - r,,| < h^. 
2 ln„ 



(18) 



Here is the smoothing length of the local evaluation point on 
the ionization front, and r^p is its radius. In this case, the tran- 
sition between the two temperature regimes is continuous (see 
Fig.|6]l. Crucially, this allows new neutral material to be ionized 
(see Fig.|7]i, and the expansion of the Hn region now follows the 
analytical solution (Eqn. |4| rather closely. 



30 



Case B 



10 



Case A 
\ 



0,1 0,3 
t (Myr) 



0,3 



Fig. 7. Growth of mass of the Hii region against time. In both 
cases we measure the mass of the gas that is completely ionized 
(ri = 10"* K). Clearly in Case B new mass is ionized while in 
Case A the mass is almost constant. 



Several other authors have encountered the same problem 
in the past. Kessel-Deynet & Burkert (2000) smooth the ioniza- 
tion fraction over a distance of the order of one local smooth- 
ing length at the position where the ionization front is located. 
Krumholz et al. (2007) considers the same problem and con- 
structs a region in which a similar procedure is followed. 
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Fig. 8. The radius of the shell is plotted against time, for both 
cases A and B along with the semi-analytical solution (Eqn. 10 1. 
The evolution of the Hn region in Case A is in substantial dis- 
agreement, while in Case B the Hn region expands according to 
Eqn. 

3.7. Computational cost 

The computational cost of locating the ionization front scales as 
Ni log (Ni), where Ni is the number of ionized SPH particles. 
This is demonstrated in Figure |9] which shows the CPU time, 
taken by the ionization subroutine, as a function of Ni, for a 
sequence of simulations performed with different total numbers 
of SPH particles, A/^q^; these are the same simulations as those 
used to test convergence in Section [4T| 

Abrupt increases in f.^^g occur when the level of refinement 
in HEALPix, {, increases by one, in response to the resolution 
criterion described in Section [J!2] Apart from this we see that 
the scaling fits Ni log (TVi) very closely, as might be expected. 

In this sequence of simulations, the overall cost of the ioniza- 
tion subroutine, over the entire duration of the simulation, scales 
as A/'spH log (N^ph), where A/^^p^ is the total number of SPH parti- 
cles. 



4. Applications 

In order to demonstrate the capabilities of the new algorithm, we 
present the results of four diff'erent applications, (i) In Section 
|4.1| we revisit the simulation described in Section [376] in which 
a spherically symmetric Hii region expands inside a spheri- 
cally symmetric, uniform-density, non-self-gravitating cloud. In 
Section [331 this simulation was used to establish the need for 
temperature smoothing. Here we use it to explore the effect of 
changing the resolution by varying the total number of SPH par- 
ticles. We show that with ^ 300,000 SPH panicles, the sim- 
ulations are well converged, both with one another, and with 
the analytic approximation derived by Spitzer (Eqn.|4]i. (ii) In 
Section |4.2| we repeat this simulation with a much more mas- 
sive cloud, and include the effect of self-gravity. The evolution 
is now well described by the semi-analytic solution derived in 
Section [2!2] (see Eqn. [TOl. (iii) In Section |4.3| we consider the 



same configuration as in Section 4. 1 except that the ionizing star 
is now initially located towards one edge of the cloud, and the 
self-gravity of the gas is neglected. As a consequence of the in- 



3x10^ particles 
6x10^ particles 
10^ particles 
2x106 particles 
NlogN 




Fig. 9. This figure shows how the algorithm scales as the number 
of ionized particles increases. The abrupt changes in f,j,yg occur 
when the level of refinement in HEALPix, {, increases by one. 
The discrepancies in t.^^^g for specific Ni are due to asymetrical 
branches of the families of rays in each simulation. 



trinsic asymmetry in the initial conditions, the Hii region bursts 
out of the cloud on one side, and the remainder of the cloud 
on the other side is accelerated away by the rocket mechanism 

we 



(Oort & Spitzer 1955), and then ionized, (iv) In Section 4.4 



consider a much lower-mass, uniform-density, self-gravitating 
cloud which is overrun by an Hii region. The ionization front 
propagating into the cloud is preceded by a shock front, which 
compresses the cloud,but does not trigger collapse. 

In all simulations, the SPH particles are initially positioned 
randomly, and then settled to produce a uniform-density "glass"; 
the temperature of the ionized gas is set to Ti = 10^ K (except 
in the boundary layer, where the temperature is smoothed to val- 
ues between these two extremes, as described in Section |3.6| i; 
the recombination coefficient into excited states only is taken to 
be ffg = 2.7 X lO"'"' cm^ s"'; and the volume occupied by the 
ionizing star is neglected. 

The parameters of all simulations are listed in Table [T| In 
this table M is the mass of the cloud, R is its radius, T„ is the 
temperature of the neutral medium, X and Y are the fractions of 
hydrogen and helium respectively, fi„ and yUj are the mean molec- 
ular weights of the neutral medium and the ionized medium re- 
spectively, //l^c is the photon rate emission of the source, D is 
the distance of the source from the centre of the cloud, R^^ is the 
initial Str0mgren radius and ATsp,, is the number of SPH particles 
used in each simulation. 



4. 1 . Spherically symmetric expansion of an Hn region in a 
non-self-gravitating uniform-density cloud 

In this application, we create a uniform-density spherical cloud 
having mass M = 1000 M©, initial radius R = 1 pc, and hence 
initial density pn = 1.6 x lO^^^gcm"-*. An ionizing source is 
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Table 1. This table is a summary of the parameters used in all applications. 



)lication 


M(Me) 


R(pc) 


r„ (K) 


X 


Y 


A'n 






£i(pc) 


^s, (Pc) 


Self-gravity 




5 3.6 




1000 


1 


10 


0.7 


0.3 


2.35 


0.678 


104« 





0.189 


OFF 


10" 


H-1 




1000 


1 


10 


0.7 


0.3 


2.35 


0.678 


1049 





0.189 


OFF 


3,6, 10,20x 10^ 


i4.2 




1.5 X 10^ 


14.6 


10 


0.7 


0.3 


2.35 


0.678 


1049 





1.43 


ON 


10^ 


H-3 




300 


1 


100 


0.7 


0.3 


2.35 


0.678 


1049 


0.4 


0.42 


OFF 


10*^ 


H.4 




20 


0.5 


100 


1 





1 


0.5 


3.2 X lO"*** 


3.5 




ON 


3x 10' 




Fig. 10. The radius of the ionization front against time, for all 
four simulations in Section 14.11 The simulations are terminated 
when the ionization front reaches the edge of the cloud. The solid 
line is the Spitzer solution, displaced upwards by AR^^- - 0. 1 pc 

to avoid confusion. Convergence is achieved for AT^py ~ 10^. 



Fig. 11. The radii of the ionization front and the shock front 



placed at the centre of the cloud, and emits ionizing photons 
at a constant rate 7Vl,c - lO'^^s"'. Using Eqn. (|2j, the initial 
Str0mgren radius is R^^ - 0.189pc. The neutral gas is assumed 
to be at Tn - 10 K. There are no gravitational forces. 

In order to establish that our code is converged, we evolve 
this configuration using different numbers of SPH particles: 



3 X 10\ 6 x 10\ 10", and 2 x 10". We terminate 



the simulations when the ionization front reaches the edge of the 
cloud. In Fig. 10 we plot the average radius of the ionization 
front, /?,p, against time for all four simulations. We see that the 
curves converge for A/'sp^ ^ 10". 

In Fig.fTTIwe plot the radius of the ionization front, R^, and 
the radius oTthe shock front, R^^, against time, for the simulation 
performed with A/^p,, =2x10^ SPH particles. We determine the 
radius of the shock front by finding the most distant particle from 
the source which has radial outward velocity > O.lkms ' 
and density p > l.lp„. The second condition is necessary be- 
cause particles near the edge of the cloud move outwards due to 
the thermal pressure gradient there, long before they are overrun 
by the shock front. The factor 1.1 is to accommodate numerical 
noise in the SPH estimate of the density of a particle. 



against time, for the simulation of Section 4.1 performed with 
A/^sp,, = 2 X 10^ SPH particles. The solid line is the semi-analytic 
solution described in Sec 12. 21 



In general, we find that the Spitzer solution (Eqn. |4|i predicts 
the radius of the ionization front well, whereas the semi-analytic 
solution (Eqn. 10 1 predicts the position of the shock front. The 
advantage of the semi-analytic solution is that it can also be used 
to treat situations in which the self-gravity of the gas is impor- 
tant, as we show in the next application. 

4.2. Spherically symmetric expansion of an Hn region in a 
self-gravitating uniform-density cloud 

In this application, we simulate a much larger, uniform-density 
spherical cloud, having mass M = 1.5 x 10^ M©, initial radius 
R - 14.6 pc, and hence initial density p„ - 7.63 x 10"^^ g cm"^. 
An ionizing source is placed at the centre of the cloud, and emits 
ionizing photons at a constant rate A/^yc - 10"*^ s"'. Using Eqn. 
Q, the initial Str0mgren radius is R^^ - 1.43 pc. The neutral 
gas is assumed to be at T„ - 10 K. We use 10'' SPH particles, 
and evolve the system for 4.5 Myr. 

In this simulation self-gravity is taken into account, but with 
the following two modifications. First, once the radius of the 
shock front has been determined (as described in |4. l| i, we ne- 
glect the gravitational acceleration of all the SPH particles out- 
side this radius. This is to prevent infall of the undisturbed neu- 
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Fig. 12. This diagram shows the evolution of the shock front 
(dashed Une) and the ionization front (dotted line) for the test 



-1.5 



described in 5 4.2 We plot Eqn. ( 10 1 for comparison (solid line) 



The evolution of shock front agrees with the prediction of our 
semi-analytical approximation. 



tral gas. Otherwise the outer parts of the cloud are already falling 
quite rapidly towards the centre by the time the shock front 
reaches them, and this complicates the interpretation of the dy- 
namics. Secondly, we only take account of the radial component 
of the gravitational acceleration. We do this to suppress frag- 
mentation of the shell, again in order to keep the dynamics of 
the shell as simple as possible (i.e. spherically symmetric). 



In Fig. 12 we plot the position of the ionization front, the 
position of the shock front, and the semi-analytic approximation. 
We see that the position of the shock front agrees very well with 
the semi-analytic approximation. 

4.3. Off-centre expansion of an Hii region 

In this application, we consider a uniform-density spherical 
cloud, having mass M - 300 Mq, initial radius R = 1 pc, 
and hence initial density pn ~ 4.85 x 10"^' gcm"^^. An ioniz- 
ing source is placed at distance D = 0.40 pc from the cen- 
tre of the cloud, and emits ionizing photons at a constant rate 
A/^yc = 10"*^ s"'; the initial Str0mgren sphere therefore has ra- 
dius Rst = 0.42 pc. The neutral gas is assumed to be at T,^ - 
100 K. The simulation uses A/^sp^ = 3x10^ particles. There are no 
gravitational forces. The simulation is terminated at f = 0.5 Myr. 



Fig. 13 shows three stages in the evolution of the cloud. At 
t ~ 0.018 Myr (Fig.[T3^), the ionization front breaks through the 
edge of the cloud on the left-hand (near) side. From this mo- 
ment on, the ionized gas can stream away freely on this side. 
However, on the right-hand side the ionization front continues 
to drive an approximately parabolic shock front into the inte- 
rior of the cloud. By f ~ 0.07 Myr (Fig. 13?), this shock has 
passed the centre of the cloud. At this stage the swept-up layer 
of neutral gas starts to break up into small "cometary knots". 
By t ~ 0.14 Myr (Fig. 13 ;), the shock front has opened up into 



a more hyperbolic shape, and has reached the right-hand edge 
of the cloud. The shocked neutral gas has spawned even more 
cometary knots. Fig. 14 shows a close-up of a few of the knots 
at t ~ 0.38 Myr. At t ~ 0.50 Myr, about 6 % of the gas remains 



3.5 



4.5 



Fig. 14. Detail of the ofF-centre expansion of the Hn region. This 
column density plot shows a fraction of the hundred cometary 
knots formed att ~ 0.38 Myr. Column density, S, is measured in 
Mq pc"^ ^= 8.7 X lO''' H2cm"^). The spatial axes {x and y) are 
labelled in parsecs. 



neutral. All this neutral gas is contained in hundreds of small 
cometary knots, which are steadily being ablated by the ionizing 
flux. 

We believe, but cannot unambiguously demonstrate, that the 
formation of these cometary knots is due to a combination of hy- 
drodynamical instabilities. In the early stages, the shell should be 
prone to the Vishniac instability (Vishniac 1983); and later on, 
when all the gas has been swept up into the shell and is being 
accelerated, it should be prone to the Rayleigh-Taylor instabil- 
ity. Therefore we expect the shell to break up into small knots. 
However, the size of the knots in our simulation is just larger 
than the resolution limit of the code, i.e. individual knots con- 
tain ~ 100 SPH particles, and therefore have masses ~ 0.1 Mq 
and diameters ~ 0.03 pc. It is not possible to establish the reality 
of the knots by performing convergence tests; if we increase the 
total number of SPH particles, Nj^, the masses of the knots de- 
crease correspondingly. We note that these knots are very rem- 
iniscent of those seen in the Helix Nebula (O'Dell & Handron 
1996), and we will attempt to model this source in a future pub- 
lication (cf. Capriotti 1973; O'Dell & Burkert 1996; Capriotti & 
Kendall 2006). 



4.4. Radiation Driven Compression 

In this application, we consider a uniform-density spherical 
cloud, having mass M - 20 Mq, initial radius R = 0.5 pc, 
and hence initial density pn ~ 2.6 x 10^^'gcm"^. We place 
an ionizing source a distance D = 3.5 pc away from the cen- 
tre of the cloud. The source emits ionizing photons at a con- 
stant rate TV. ^ - 3.2 x 10"^^ s"'. Hence the number-flux of ion- 

LyC 

izing photons incident on the near side of the cloud is O ~ 
2.97 X 10^ cm"^ s Since D '» R, the rays are parallel and the 
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Fig. 13. Column-density plots of the off-centre expansion of the Hii region. Column-densities, S are measured in 
Mopc"^ 8.7 X lO'^Hacm"^). The spatial axes {x and y) are labelled in parsecs, and times in Myr are given in the top right- 
hand corner of each frame. 



flux is constant. In this one application we assume X - \ instead 
of X = 0.7, i.e. pure hydrogen with yUn = 1 and yUi = 0.5. The 
neutral gas is assumed to be at - 100 K. The sound speeds 
are therefore Cn - 0.9 km s"^ and Ci = 12.8 km s"'. The simula- 
tion uses A^sp,, = 3 X 10^ particles, and self-gravity is taken into 
account. The simulation is terminated at f = 2.5 Myr. 

The above initial conditions are chosen to reproduce - as 
closely as possible - the initial conditions of the 2D model by 
Lefloch & Lazareff" (1994). We note that from the outset the 
cloud is over-pressured, and left to its own devices would simply 
disperse on a timescale of ~ 1 Myr. 

Fig. 15 shows the evolution of the cloud. The ionizing flux 



propagates upwards and rapidly boils off the outer layers on the 
near side of the cloud. At f ~ 0.036 Myr (Fig. 15 i) a shock front 
starts to compress the remaining neutral gas. At the same time, 
the north hemisphere starts to expand due to the thermal pres- 
sure of the atomic hydrogen. At f ~ 0.13 Myr (Fig. 15 i) a dense, 
prolate, approximately ellipsoidal core forms. By f ~ 0.18 Myr 
(Fig. 15 ;), the prolate core has semi-major axis 0.08 pc and 
semi-minor axis 0.02 pc; its mass is ~ 6M0. It is therefore ther- 
mally sub-critical, and does not collapse to form a star. Instead 



151) 



By 



it is steadily ablated by ionization. By f ~ 0.21 Myr (Fig 
the remnants of the cloud start to develop a cometary tail 
t ~ 2.5 Myr (not shown on Fig. 15 1, the last vestiges of the cloud 
are ionized, and they are ~ 28 pc from the ionizing star. 

Fig.[T6]plots the total mass of neutral gas against time. The 
undulations seen at f ~ 0.25 Myr, and at f ~ 0.5 Myr are acous- 
tic oscillations, excited as the cloud responds to the increase in 
external pressure. 

Gritschneder et al. (2009; hereafter G09) have also simulated 
radiatively driven compression. However, the cloud that G09 
treat (with mass 96 M©, radius 1.6 pc, and temperature 10 K) 
is much more massive and much colder than the one we treat 
here (20 Mq, 0.5 pc, 100 K), and therefore it is much less resis- 
tant to compression and more prone to triggered gravitational 



collapse. Furthermore, in their first two simulations, G09 use a 
larger ionizing flux than we do. As a consequence, the clouds in 
their simulations are more strongly compressed than ours, par- 
ticularly at the leading edge (i.e. the edge exposed directly to 
the ionizing flux), and are triggered into gravitational collapse. 
In contrast, our cloud is more mildly compressed, and evolves 
towards a centrally condensed configuration but never becomes 
gravitationally unstable. 



5. Discussion & Conclusions 

We have introduced a new technique for treating the propaga- 
tion of ionizing radiation in SPH simulations of self-gravitating 
gas dynamics. The method uses the HEALPix algorithm to tes- 
sellate the celestial sphere, and solves the equation of ionization 
equilibrium along the rays associated with each tessera; rays are 
split hierarchically to produce greater resolution wherever it is 
required, i.e. to ensure that the resolution of the radiation trans- 
fer matches the resolution of the hydrodynamics, locally. This 
makes the algorithm very computationally efficient. 

The algorithm has been incorporated into the new Smoothed 
Particle Hydrodynamics code SEREN (Hubber et al., in prepa- 
ration), and tested against a number of known analytic and semi- 
analytic problems. These tests are presented here. The code fol- 
lows the expansion of a spherically symmetric Hii region in a 
non-self-gravitating gas (Spitzer solution); the expansion of a 
spherically symmetric Hii region in a self-gravitating gas (new 



semi-analytic solution described in Sec 2.2 1; the rocket accelera 



tion and subsequent ablation of a massive cloud irradiated by an 
ionizing star on one side (Oort & Spitzer 1955); and the radia- 
tively driven compression of a pre-existing dense core engulfed 
by an expanding Hii region. The code will be used in future to ex- 
plore the role of Hii regions in triggering and regulating star for- 
mation, injecting turbulent energy into the interstellar medium, 
and eroding molecular clouds. 
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Fig. 15. Column-density plots of radiation driven compression. Column-densities, S are measured in Mq pc"^ (= 1 .3 x 10^*^ cm"^). 
The spatial axes (x and y) are labelled in parsecs, and times in Myr are given in the top right-hand corner of each frame. 




I (Myr) 



Fig. 16. This diagram shows the evaporation in time of the neu- 
tral gas. The solid line represents the simulation discussed in 

5 4.4 whereas the dash line the simulation by Lefloch & Lazareff 
(1994). The discrepancy between our simulation and the simu- 
lation by Lefloch & Lazareff" (1994) at time t > 1.1 Myr is due 
to the movement of the cloud further from the star, where the 
radiation flux drops in our model but stays constant in Lefloch 

6 Lazareff (1994) plane parallel model. 
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Appendix A: Resolution requirements to resolve 
the ionization front 

Consider an Hn region having uniform density p\. The SPH par- 
ticles have a universal mass m^p^. If each SPH particle has N^^^^ 
neighbours, then the diameter, d, of an SPH particle (i.e. the di- 
ameter of its smoothing kernel) is given by 



n<f Pi 



6M„ 



1/3 



(A.l) 



Since d is in effect the resolution of the simulation, the ioniza- 
tion front can only be resolved \f d % A/?,p. Using Eqn. (|3]l, this 
requirement reduces to 



(A.2) 



Combining Eqns. (|6]l and (|5]), we see that the mass of the Hn 
region is given by 



Mi = 



«B Pi 



(A.3) 



The number of SPH particles required to model the Hn region is 
therefore 



Mi 

6 0-^A^N.E,BM,yC Pi 

(20)-' 7T m 



^ 5 X 10" 



50 



1049 s-1 



10-2"gcm-3 



(A.4) 
(A.5) 
(A.6) 



Evidently this is a prohibitive requirement on Ni, all the more 
so when one allows that the mass of neutral gas is likely to be 
even greater than the mass of ionized gas, and so this will require 
additional SPH particles. 



